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Collinear Earth-Moon libration points have emerged as locations with immediate 
applications. These libration point orbits are inherently unstable and must be maintained 
regularly which constrains operations and maneuver locations. Stationkeeping is 
challenging due to relatively short time scales for divergence, effects of large orbital 
eccentricity of the secondary body, and third-body perturbations. Using the Acceleration 
Reconnection and Turbulence and Electrodynamics of the Moon’s Interaction with the 
Sun (ARTEMIS) mission orbit as a platform, the fundamental behavior of the trajectories 
is explored using Poincare maps in the circular restricted three-body problem. Opera- 
tional stationkeeping results obtained using the Optimal Continuation Strategy are 
presented and compared to orbit stability information generated from mode analysis 
based in dynamical systems theory. 
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1. Introduction 

Earth-Moon collinear libration points have emerged as 
locations with immediate applications. The selection of 
Earth-Moon (EM) libration point orbit orientations and 
amplitudes, as well as the inherent stationkeeping chal- 
lenges, are correlated and achieving the best result is 
nontrivial. Thus, this paper offers the Acceleration Recon- 
nection and Turbulence and Electrodynamics of the Moon’s 
Interaction with the Sun (ARTEMIS) mission as a platform 
to explore: the theory of EM libration point orbit evolution, 
the required modeling of the libration point orbits, an 
assessment of Earth-Moon libration point orbit stability. 


* This paper was presented during the 1st DyCoSS in Porto. 

* Corresponding author at: 8800 Greenbelt Rd., Greenbelt, MD 20771- 
2400, USA. Tel.: +1 301 286 6082; fax: +1 765 494 0307. 

E-mail addresses: david.cfolta@nasa.gov (D.C. Folta), 
tpavlak@purdue.edu (T.A. Pavlak), ahaapala@purdue.edu (A.F. Haapala), 
howeII@purdue.edu (ICC Howell), 
mark.a.woodard@nasa.gov (M.A Woodard). 


the use of Poincare maps to define and, potentially, select 
orbit characteristics, and in-flight stationkeeping strate- 
gies. Results using operational data are demonstrated in 
orbit analysis using Poincare maps, and in the implemen- 
tation of stationkeeping strategies as applied to ARTEMIS. 

The libration point orbits employed for ARTEMIS are 
inherently unstable and stationkeeping maneuvers applied 
approximately once per week are required. The relatively 
short time scales for divergence, due to the effects of large 
orbital eccentricity of the secondary and solar gravitational 
and radiation pressure perturbations, prove challenging for 
orbital stationkeeping. However, the Optimal Continuation 
Strategy (OCS), which employs various numerical methods 
to determine maneuver locations and optimize stationkeep- 
ing delta-V (AV), successfully minimizes fuel cost during 
ARTEMIS mission operations while offering quality naviga- 
tion tracking and maneuver planning scenarios. To examine 
libration point orbit control methods, the operationally 
demonstrated and validated Optimal Continuation Strategy 
is compared with an implementation of Floquet mode 
information as calculated from navigation states. The OCS 
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yields the optimal maneuver locations and AV directions for 
the given mission constraints and it is demonstrated that 
these AV directions compare closely with the stable eigen- 
vector directions computed via Floquet analysis. 

2. Background 

2.1. The ARTEMIS mission 

As the first Earth-Moon libration point orbiting mis- 
sion, ARTEMIS serves as the basis (investigative platform) 
for this exploration. More specifically, ARTEMIS is the first 
mission flown to and continuously maintained in orbit 
about both collinear Earth-Moon libration points, EM Li 
and EM l 2 11—5]. Originally, five spacecraft comprised the 
successful Time History of Events and Macroscale Inter- 
actions during Substorms (THEMIS) mission. The ARTE- 
MIS mission was characterized by the transfer of the two 
outermost THEMIS spacecraft from their elliptical Earth 
orbits and, with lunar gravity assists, re-directed them 
to the vicinity of EM L 2 via transfer trajectories that 
exploited the Sun-Earth multi-body dynamical environ- 
ment. Two identical ARTEMIS spacecraft, named PI 
and P2, entered the space near the Earth-Moon libration 
point orbits in 2010 on August 25 and October 22, 
respectively. Once the Earth-Moon libration point orbits 
were achieved, they were maintained for approximately 


eleven months. The PI spacecraft was inserted into orbit 
in the vicinity of EM L 2 and, during the stationkeeping 
phase, was transferred from EM L 2 to EM Li. The P2 
spacecraft, however, inserted near EM L 2 , but immediately 
transferred into an EM Li orbit where it was maintained 
for the remainder of the mission. At the conclusion of the 
stationkeeping phase, both spacecraft departed their EM 
Li orbits and inserted into elliptical, low inclination lunar 
orbits in 2011 on June 27 and July 17, respectively. 

Various views of the EM L 2 and Li libration point orbits 
associated with the ARTEMIS PI spacecraft appear in 
Figs. 1 and 2. Similarly, projections of the ARTEMIS P2 
EM Li libration point orbit are depicted in Figs. 3 and 4. 

There were no size or orientation requirements on 
these orbits other than to minimize the libration point 
orbit insertion and stationkeeping AV costs and, ulti- 
mately, to allow for transfer into low inclination lunar 
orbits. Both ARTEMIS spacecraft had limited (combined) 
deterministic and statistical stationkeeping AV budgets of 
15 m/s and 12 m/s for PI and P2, respectively. These AV 
budgets included any libration point orbit stationkeeping 
requirements as well as the transfers between libration 
points and into lunar orbit. The amplitudes of the ARTE- 
MIS orbits were governed by the EM L 2 orbit insertion 
conditions as dictated by the ballistic Sun-Earth to Earth- 
Moon transfers that were utilized. The resulting PI and P2 
EM Li orbit y-amplitudes were approximately 60,000 km 





Fig. 2. ARTEMIS PI libration point orbits, additional projections, (a) x-z projection and (b) y-z projection. 
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with the PI L 2 orbit y-amplitude near 68,000 km. 
A complete list of the general dimensions for the evolving 
ARTEMIS libration point orbits appears in Table 1 . 


2.2. Environmental modeling 

The ARTEMIS PI and P2 trajectories, like many colli- 
near libration point orbits, are dynamically unstable and 
even small perturbations can significantly affect their 
evolution. During operations, if any inaccuracies in the 
dynamical model are present, a spacecraft can depart 
from the desired EM L, or L 2 orbit along an unstable 
manifcld, either towards the Moon or along an escape 
trajectory towards the Earth or, even, away from the 
Earth-Moon system. The perturbations required to affect 
these changes are so small that even mis-modeled accel- 
erations due to solar radiation pressure or additional 
gravitating bodies can result in departure from the nom- 
inal libration point orbit. Thus, for successful mission 
operations, the system must be modeled as a true four- 
body problem, including the Sun’s gravity as an important 
third-body perturbation; accelerations due to lunar 
eccentricity and solar radiation pressure are incorporated 
as well. The relative effect of these force terms on the 
evolution of the ARTEMIS P2 trajectory are explored by 
Pavlak and Howell [6]. Operationally, the ARTEMIS mis- 
sion trajectory design and analysis utilized a full DE421 
Moon-Earth-Sun ephemeris model that also incorporated 
solar radiation pressure acceleration based upon the 
spacecraft mass and cross-sectional area - in this case, a 
simplified cannon ball model was used - and a gravity 
potential model for the Earth with degree and order eight. 
The high-fidelity orbit propagation was conducted using 
both Goddard Space Flight Center’s General Mission 
Analysis Tool (GMAT) and AGI’s STK/Astrogator software 
package. Successful libration point orbit stationkeeping 
operations rely on high-fidelity modeling to accurately 
predict trajectory evolution and, ultimately, to design AV 
maneuvers that ensure the mission objectives are satisfied. 


Table 1 

ARTEMIS libration point orbits parameters. 


Orbit parameter 

PI at l, 

PI at L 2 

P2 at L, 

Maximum x-amplitude (km) 

23,656 

32.686 

30.742 

Maximum y-amplitude (km) 

58,816 

63.520 

67,710 

Maximum z-amplitude (km) 

2387 

35198 

4680 

Orbital period (days) 

13.51 

15.47 

14.19 


3. Theoretical background and libration point orbit 
analysis 

While high-fidelity modeling is important for final 
mission design and operations, the circular restricted three- 
body (CR3B) problem is a simplified multi-body model 
that can provide much insight for pre-flight trajectory 
planning and offers a very useful perspective during 
post-mission analysis for ARTEMIS. The libration point 
orbits for the ARTEMIS mission were selected to minimize 
libration point orbit insertion and stationkeeping costs 
and to enable low-cost transfers to low inclination lunar 
orbits; there were no specific size or orientation require- 
ments on these orbits. Analysis of the libration point 
orbits from the perspective of the CR3B problem via 
Poincare maps offers further insight into the orbital 
evolution of the ARTEMIS spacecraft trajectories and 
dynamical evidence that demonstrates why the selected 
trajectories were ultimately successful in satisfying the 
mission constraints. 

3.1. The circular restricted three-body problem 

In the CR3B model [7], the motion of a spacecraft, 
assumed massless, is governed by two gravitationally 
massive primary bodies, the Earth and the Moon, each 
represented as a point mass. The orbits of the primary 
bodies are assumed circular relative to the system bar- 
ycenter. A baiycentered rotating frame is defined such 
that the rotating x-axis is directed from the Earth to the 
Moon, the z-axis is parallel to the direction of the angular 
velocity of the primary system, and the y-axis completes 
the right-handed triad. The spacecraft state relative to the 
Earth-Moon barycenter is defined in terms of rotating 
coordinates as x=[x,y,z,x,y,i]. Note that bold symbols 
denote vector quantities. The mass parameter, p, is 
defined 


where m, and m 2 correspond to the mass of the Earth 
and Moon, respectively. The first-order, nondimensional, 
vector equation of motion is 

*=/(*) ( 2 ) 

where the vector field, /(x), is defined 

f(x) = [x,y,z,2y+U x .-2x + U y ,U z ] (3) 

The pseudo-potential, U, is defined U(x,y,z) = ^ + £ + 
£(x 2 +y 2 ) with the nondimensional Earth-spacecraft 
and Moon-spacecraft distances written as di and d 2 , 
respectively. The quantities U x ,U y ,U z represent partial 
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derivatives of U with respect to rotating position coordi- 
nates. The single, scalar integral of the motion, known as 
the Jacobi constant, C, is represented as C = 2 U-v 2 , where 
v = (x 2 +y 2 +z 2 ) 1/2 . 

Five equilibrium points exist, including three collinear 
libration points, L U L 2 , and I 3 , that lie along the x-axis, and 
two equilateral points U and L5. Linear analysis of the 
collinear points [7-10] reveals that they possess a topo- 
logical structure of the type saddle x center x center. 
Thus, asymptotic flow to and away from the libration 
points is possible via the stable and unstable manifolds, 
respectively; periodic and quasi-periodic orbits exist 
within the center subspace. 

3.2. Libration point orbits and stability 

A linear analysis of the flow in the vicinity of the 
libration points provides insight into the evolution of 
libration point orbits in the nonlinear CR3B problem. 
Linearization of the equations of motion about the colli- 
near libration points allows for states to be selected such 
that they exist within the center subspace, yielding the 
following variational equations describing orbits relative 
to the libration point: 

<5x(t)=A< cos(yt+<f>), (4) 

Sy(t) = A y sin(yt+0), (5) 

Sz(t) = A z sin(cot -h/0, ( 6) 

where iy. ico represent the complex eigenvalues asso- 
ciated with the local planar and out-of-plane center 
manifolds, respectively; angles <f>, \j/ represent phase 
angles and amplitudes A x , A y are related by a proportion- 
ality constant. Because y^eo, these equations describe 
quasi-periodic motion in the vicinity of the collinear 
points whenever A x ^0 and A Z A 0. Selecting A z —0 yields 
the planar, periodic Lyapunov orbits while choosing 
A x = Ay — 0 produces the periodic vertical orbits. An 
example of a quasi-periodic orbit in the linear model 


appears in black in Fig. 5(a); note that the orbit is 
linearized relative to the Li point, but is plotted in a 
Moon-centered view. The z-amplitude evolution appears 
in Fig. 5(b), and illustrates the constant amplitude, A z . 

By adjusting the phase angles <f>, i/c it is possible to 
enter the libration point orbit in different locations [11], 
As an example, the red segments in Fig. 5 represent shifts 
in angles <f> and \j/ that allow for entry into the orbit at a 
desired point within the z-amplitude evolution cycle. 

The variational equations exist within the framework of 
a linear analysis, however, periodic and quasi-periodic orbits 
also exist in the full nonlinear model, described by Eqs. (2) 
and (3) [12], Several methods exist to compute periodic and 
quasi-periodic libration point orbits with the inclusion of 
nonlinear effects [13-17J. A quasi-periodic orbit in the fully 
nonlinear CR3B problem is depicted in Fig. 6(a); the corre- 
sponding z-amplitude evolution appears in Fig. 6(b). 

Clearly, the z-amplitude is no longer constant, but cycles 
between high and low z-amplitude modes. By shifting the 
entry phasing in the nonlinear model, it is possible to enter 
the quasi-periodic orbit at a different location and alter the 
z-amplitude mode. The red highlighted region of Fig. 6 
illustrates how the entry location in the orbit can be shifted 
to enter during the nearly planar mode. The entry clearly 
influences the z-amplitude evolution later along the trajec- 
tory and/or the maneuver history required to maintain 
the orbit 

To determine the orbit stability properties of a periodic 
libration point orbit the eigenstructure, i.e., eigenvectors 
and eigenvalues, is analyzed. The monodromy matrix, M, is 
generated by integrating the state transition matrix (STM), 
<£(f2,ti), along with the trajectory state vector for one 
orbital period, T, i.e., M = 4>(t+T,t). By computing the 6 
eigenvalues, of the monodromy matrix M, one can 

determine the stable eigenvalue^) (|2j| < 1), unstable 
eigenvalue(s) (]A,| > 1), and eigenvalues reflecting behavior 
within the center subspace (J2 ( J=1) associated with a 
particular orbit. An orbit for which all |A ( | = 1 is considered 
marginally stable; those orbits for which at least one pair of 
stable and unstable eigenvalues exist are termed unstable. 


a b 



Time (Days) 


Fig. S. Quasi-periodic orbit.and corresponding z-amplitude evolution in the linear model, (a) Quasi-periodic orbit in linear model. Moon-centered view 
and (b) z-araplitude evolution over time. 
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Fig. 6. Quasi-periodic orbit in the CR3B model cycles through high and low z-amplitude modes, (a) Qtiasi-periodic orbit in CR3B model, Moon-centered 
view and (b) z-amplitude evolution over time. 
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Fig. 7. Poincare map depicting periodic and quasi-periodic libration point orbit structures in the vicinity of in the Earth-Moon system for C=3.08O. 
(a) Poincare map corresponding to crossings of the x-y plane, with sample orbits featured and (b) z-amplitude evolution of selected orbits. 


3.3. Poincare maps 

For a more complete picture of the available libration 
point orbits at a particular energy level in the CR3B 
problem, and for potential design options, it is useful to 
employ Poincare maps. Through the use of a Poincare 
map, an n-dimensional continuous-time system is reduced 
to a discrete-time system of (n-l)-dimensions. By addi- 
tionally constraining the Jacobi constant, C, the problem is 
reduced to (n-2)-dimensions and, thus, the map for the 
CR3B problem is represented in 4-D. To generate a planar 
Poincare map, a surface-of-section, S, is defined such that 
I is transverse to the flow. A commonly used I is one that 
represents a surface-of-section corresponding to crossings 
of the x-y plane. To compute the map. trajectories are 
integrated using Eqs. (2) and (3), and crossings of I are 
recorded and displayed. 

Consider the map in Fig. 7 as projected into the x-y 
plane, produced to resemble the maps demonstrated by 


Gomez et al. [18], as well as Kolemen et al. [16] For the 
selected value of C, several periodic orbits exist, including 
a planar Lyapunov orbit (green), a vertical orbit (dark 
blue), and the northern and southern halo orbits. In the 
x-y projection, the halo orbits share the same crossings of 
the map; the northern halo is featured in magenta in 
Fig. 7(a). Utilizing multiple shooting and 30-50 “stacked 
revolutions” of the periodic orbits, quasi-periodic regions 
of the Poincare map are populated using a numerical 
continuation scheme that begins at the map crossing 
associated with a periodic orbit and steps outward to 
yield successively larger quasi-periodic trajectories. Sur- 
rounding the periodic vertical orbits are quasi-periodic 
orbits, often denoted Lissajous orbits, that exist within 
the center subspace associated with the vertical orbit. A 
sample Lissajous trajectory is featured in cyan. Similarly, 
the quasi-halo orbits lie in the center manifold of the 
central halo orbit. Examples of small and large northern 
quasi-halo orbits appear in red and orange, respectively. 
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The z-amplitude evolutions corresponding to the large 
quasi-haio and Lissajous orbits featured in Fig. 7(a) appear 
in Fig 7(b). The periodic halo and vertical orbits possess 
constant amplitudes A z , whereas the quasi-periodic orbits 
exhibit oscillating values of A z . The crossings of the Lissajous 
orbits are contained within the central region of the map; 
therefore, these orbits do not possess the nearly planar 
modes demonstrated in Fig. 6 that are desired for ARTEMIS- 
type libration point trajectories. The quasi-halo orbit cross- 
ings occur in the upper and lower regions of the map, thus, 
both high and low z-amplitude modes are facilitated by 
selection of a quasi-halo orbit for this value of C. These 
distinct regions of quasi-periodic behavior are explored in 
detail by Barden and Howell [19]. 

3.4. ARTEMIS Poincare maps 

To gain insight into the types of orbits selected for 
ARTEMIS, Poincare maps are generated to display the 
specific orbit structures associated with the libration point 
orbit energy levels associated with the PI and P2 spacecraft 
while in libration point orbits. The ARTEMIS trajectories 
were designed with higher-fidelity ephemeris modeling, 
and the true paths possess discontinuities in the form of 
small AV maneuvers. Thus, to analyze the libration point 
orbits using maps, it is desirable to compute orbits in the 
CR3B model that are qualitatively similar to those in the 
ARTEMIS mission, as depicted in Figs. 1-4. Analogues to the 
ARTEMIS PI and P2 trajectories are computed in the CR3B 
problem as explained by Pavlak and Howell [20]. To 
summarize the approach, an initial guess is developed by 
first “sacking” planar Lyapunov orbits for a desired number 
of revolutions about L 2 and/or L t . Stable and unstable 
manifold segments are incorporated to transfer between 
L 2 anc L- libration point orbits and, ultimately, from the 
vicinity of L, to the Moon. A multiple shooting algorithm is 
then utilized to compute continuous, end-to-end trajec- 
tories in the CR3B problem with starting points that are 
consistent with the ARTMEMiS L 2 orbit insertion states 


obtained from operational navigation data. For each of the 
converged trajectories, the Jacobi value, C, is evaluated and 
maps are generated that correspond to the appropriate 
libration point and the specified energy level. Because the 
ARTEMIS orbits and the converged CR3B libration point 
orbits appear qualitatively the same, the results from the 
discussion of the CR3B orbits can be applied to the actual 
ARTEMIS orbits as well. 

Maps for the PI L 2 , PI Li, and P2 Li libration point orbits 
appear in Figs. 8-10, respectively. Clearly, each of the three 
ARTEMIS libration point orbits possesses x-y plane crossings 
(in red) that lie in the outermost quasi-halo region of the 
map, indicating that each of the ARTEMIS orbits is a south- 
ern large quasi-halo orbit To demonstrate the long-term 
evolution of the ARTEMIS trajectories, quasi-halo orbits with 
map crossings that lie close to each of the ARTEMIS libration 
point orbit crossings are selected and appear in Figs. 8(b), 
9(b), and 10(b), in addition to the crossings associated with 
the converged ARTEMIS CR3B trajectories plotted in red. The 
selected quasi-halo orbit crossings are highlighted in color 
on the maps, and clearly lie close to the ARTEMIS spacecraft 
crossings. For both the ARTEMIS PI and P2 trajectories, 
recall that the libration point orbit insertion conditions are 
dictated by the ballistic Sun-Earth to Earth-Moon transfers 
that were utilized and that the orbits were selected, in part, 
to minimize the AV required to transfer between L 2 and Li 
and, eventually, into low inclination lunar orbits. Analysis of 
the large quasi-halo orbits demonstrates that these orbits 
supply both high and low z-amplitude phases. Note that the 
Lissajous family of quasi-periodic orbits does not provide 
solutions with the near-planar phases that are required 
based on the ARTEMIS mission constraints. Thus, the 
ARTEMIS P2 L 2 trajectory is consistent with the large 
quasi-halo trajectories because they facilitate both a highly 
out-of-plane insertion condition and a nearly planar depar- 
ture to the EM Li side of the Moon. After departing the L 2 
quasi-halo within the nearly planar phase, PI enters a large 
Li quasi-haio, but employs only the low A* phase of the orbit 
to allow for the eventual transfer into a low inclination 


a b 



Fig. 8. Poincare map associated with the ARTEMIS PI t 2 orbit (C=3.105). (a) Map corresponding to PI i 2 orbit and (b) PI f 2 orbit (red) with southern 
quasi-haio of similar size (purple). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this 
article.) 
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a b 
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Fig. 9. Poincare map associated with the ARTEMIS PI L , orbit (C=3.105). (a) Map corresponding to PI L, orbit and (b) PI L, orbit (red) with southern quasi- 
halo of similar size (green). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.) 


a b 



Fig. 10. Poincare map associated with the ARTEMIS P2 L, orbit (C= 3.080). (a) Map corresponding to P2 L, orbit and (b) P2 L, orbit (red) with southern quasi- 
halo of similar size (blue). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.) 


lunar orbit. The arrival conditions provided by the ballistic 
Sun-Earth to Earth-Moon transfer phase for the ARTEMIS P2 
spacecraft possess a much smaller out-of-plane component 
when compared with PI, however, the requirement to 
transfer to a low inclination lunar orbit remains. Thus, 
ARTEMIS P2 also exploits the low A r phase of the large L 7 
quasi-halo to reduce the AV costs associated with insertion 
into a low inclination lunar orbit. 

4. Libration point orbit stationkeeping strategies 

A variety of stationkeeping strategies have previously 
been investigated for applications in the Sun-Earth system 
and near the Earth-Moon libration points. To be operation- 
ally useful for the ARTEMIS mission, a stationkeeping 
strategy must satisfy several conditions: utilize high- 
fidelity ephemeris models, yield optimal solutions, and be 


applicable to Earth-Moon Li or L 2 libration point orbits and 
any transfer between them. There were no specific con- 
straints on the size or evolution of the ARTEMIS trajectories 
as long as the mission constraints were satisfied so station- 
keeping strategies that require strict adherence to a refer- 
ence trajectory were not well-suited to this applicatioa 
Ultimately, the stationkeeping approach selected for the 
ARTEMIS mission was able to incorporate high-fidelity 
modeling of the Earth-Moon dynamical environment and 
accommodate spacecraft constraints while meeting mission 
objectives. 

4. 1 . Previous stationkeeping strategies 

Numerous references in the literature offer discussions 
of stability and control for vehicles at both collinear and 
triangular libration point locations. Hoffman [21] and 
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Farquhar [22] both generate analysis and discussions for 
stability and control of spacecraft in Earth-Moon collinear 
Ei and l 2 locations, respectively, within the context of 
classical control theory and/or linear approximations. 
Renault and Scheeres offer a statistical analysis approach 
[23]. Howell and Keeter [24] address the use of Floquet 
theory to select stationkeeping maneuvers to eliminate 
the unstable modes associated with a reference orbit. 
Gomez et al. [25] independently develop and apply the 
Floquet mode approach specifically to translunar libration 
point orbits. Marchand and Howell [26] discuss stability 
including the eigenstructures near the Sun-Earth loca- 
tions Folta and Vaughn [27] present an analysis of 
staticnkeeping options and transfers between the Earth- 
Moon locations, and the use of numerical models that 
include discrete linear quadratic regulators and differen- 
tial correctors. Pavlak and Howell [20] demonstrate a 
long-term orbit maintenance technique incorporating 
multiple shooting. Lastly, Folta et al. [1,2,4], Woodard 
et al. [3], and Sibeck et al. [5] provide both a review of 
various Earth-Moon libration point orbit stationkeeping 
methods, as well as detailed operational stationkeeping 
and transfer results for the ARTEMIS mission. 

4.2. Optimal Continuation Strategy (OCS) 

Based on a comparison of existing libration point orbit 
stationkeeping methods by Folta et al. [1,2] and the 
imposed operational constraints, it was determined that 
the Optimal Continuation Strategy (OCS) offered a flexible 
method for maintaining the ARTEMIS Earth-Moon libra- 
tion point orbits. Fundamentally, the OCS is designed to 
maintain the spacecraft in the vicinity of a libration point 
orbit for 1-2 revolutions downstream. The OCS achieves 
this objective by computing maneuvers such that the 
spacecraft satisfies a set of user-defined constraints at 
successive crossings of the x-z plane. Often, these con- 
straints are expressed in terms of either x-position or 
x-velccity. At its foundation, the OCS is formulated as a 
two-point boundary value problem and can be solved 
using any appropriate differential corrections scheme. 
Nonlinear optimization such as sequential quadratic pro- 
gramming (SQP) can also be incorporated to yield locally 
optimal stationkeeping maneuvers. Initially, the algo- 
rithm computes a maneuver to achieve constraints one 
x-z plane crossing downstream. Once a converged solu- 
tion is obtained, the targeting procedure is repeated to 
satisfy constraints at successive future x-z plane crossings 
- typically 3-4 crossings maximum. The OCS is imple- 
mented in this manner at all subsequent stationkeeipng 
maneuver locations to maintain a libration point orbit for 
an arbitraiy number of revolutions without requiring a 
reference solution. 

5. Stationkeeping operations for the ARTEMIS mission 

5.2. ARTEMIS OCS implementation 

The flexible nature of the Optimal Continuation Strat- 
egy and the fact that it does not require a reference 
solution, made it well-suited to ARTEMIS stationkeeping 


operations. In this implementation of the OCS, constraints 
on x-velocity were incorporated at x-z plane crossings. 
However, the specific x-velocity targets employed in this 
application differed slightly between the EM L 2 and Li 
orbits and were determined empirically during pre- 
mission analysis [2]. The velocity targets to maintain the 
PI spacecraft in orbit about EM L 2 consisted of two 
different x-velocity targets, depending on the spacecraft 
location in the orbit, i.e., the side closer or farther from the 
Moon. For x-z plane crossings on the far side of the L 2 
orbit (farther from the Moon), the OCS targeted an 
x-velocity -20 m/s with a tolerance of 1 cm/s. Crossing 
targets on the near side of the L 2 trajectory (nearer to the 
Moon), however, required an x-velocity +10 m/s with a 
tolerance of 1 cm/s. For both the PI and P2 orbits in the 
vicinity of EM Et, the x-velocity, x, at each x-z plane crossing 
was constrained such that -10 cm/s < x < 10 cm/s. High- 
fidelity orbit propagation incorporating the dynamical mod- 
els discussed in Section 2.2 was performed using explicit 
Runge-Kutta-Vemer 8(9) and a Dormand-Prince 8(9) inte- 
grators in GMAT and STK/Astrogator, respectively. Each 
stationkeeping Ay maneuver was minimized using SQP 
optimizers including VF13AD from the Harwell library and 
Matlab’s fmincon. To ensure that the spacecraft remained in 
the vicinity of the Earth-Moon libration point orbits in the 
near-term, up to four downstream crossings were used 
during the OCS targeting procedure and stationkeeping 
maneuvers were executed twice per orbit - approximately 
once per week. 

5 2. Operational ARTEMIS stationkeeping results 

The OCS was tailored to the ARTEMIS Earth-Moon 
libration point orbits and was used to successfully main- 
tain the PI and P2 spacecraft in the vicinity of EM L, and 
L 2 during the 11 months of mission operations. The 
operational stationkeeping maneuvers performed for the 
ARTEMIS PI and P2 spacecraft are summarized in Table 2. 
The actual chronological Ay histories for the stationkeep- 
ing maneuvers executed for both the ARTEMIS PI and P2 
spacecraft are depicted in Fig. 11. The AV for each 
maneuver for the PI spacecraft appears in Fig. 11(a), with 
the estimated annual maintenance cost for PI in E 2 and E, 
orbit presented in Fig. 11(b) and (c), respectively. Note 
that, for each maneuver, the estimated annual mainte- 
nance cost is extrapolated for one year based on the 
operational ARTEMIS maneuver schedule and the total 
stationkeeping Ay cost accumulated to date. The station- 
keeping Ay history for ARTEMIS P2 in the EM Ei orbit is 
depicted in Fig. 11(d). The estimated annual stationkeep- 
ing costs for P2 before and after a change in the coefficient 


Table 2 

ARTEMIS PI and P2 stationkeeping results. 


AV cost 

PI at L 2 

PI at Li 

P2 at Li 

Total AV (cm/s) 

244.0 

155.0 

324.0 

Minimum AV (cm/s) 

6.96 

1.17 

133 

Maximum AV (cm/s) 

22.64 

27.90 

37.89 

Mean AV (cm/s) 

13.51 

7.21 

10.85 

Standard deviation (cm/s) 

5.44 

7.60 

10.31 
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Fig. 11. Operational stationkeeping costs for the ARTEMIS mission, (a) PI individual SK AV, (b) Pt EM i 2 libration point orbit estimated annual AV, 
(c) PI EM 1. 1 libration point orbit estimated annual AV, (d) P2 individual SK AV maneuver, (e) P2 EM £., libration point orbit estimated annual AV (pre-C r 
change) and (f) P2 EM L, libration point orbit estimated annual AV (post-Q change). 


of reflectivity, C r , appear in Fig. 11(e) and (f), respectively. 
Initially, a constant value of C n obtained from analysis 
conducted before libration point orbit insertion, was utilized 
for ARTEMIS P2 stationkeeping activities. However, from 
stationkeeping maneuver #22 onward, an updated value of 
C r based on libration point orbit navigation data was 
implemented; the predictive capability of the dynamical 
models was improved and, consequently, the stationkeeping 
costs for the remainder of the ARTEMIS P2 trajectory were 
significantly reduced. For both the PI and P2 spacecraft, the 
general decrease in the estimated annual maintenance cost 
is attributed to the improved modeling and onboard execu- 
tion of stationkeeping maneuver thrust arcs. 

A summary of stationkeeping costs for ARTEMIS PI 
and P2 following libration point orbit insertion, not 
including the axial correction maneuvers to extend the 
mission an additional three months [1], are listed in 
Table 3. Recall that the AV costs per year are extrapolated 
based on the operational ARTEMIS maneuver schedule 


Table 3 

ARTEMIS stationkeeping summary. 


AV cost 

ARTEMIS PI 

ARTEMIS P2 

Total AV (m/s) 

3.99 

3.24 

Projected annual AV, U (m/s) 

7.39 

N/A 

Projected annual AV, L t (m/s) 

5.28 

5.09 


and costs observed during the eleven month mission. 
More detailed ARTEMIS operational stationkeeping AV 
data, including maneuver epoch information, is provided 
by Folta et al. [1]. 

5.3. ARTEMIS post-mission mode analysis 

From the theoretical perspective, libration point orbit 
stationkeeping research has been ongoing for over a 
decade [23-25,28]. One particular investigation 
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applicable to this analysis is a collaborative effort 
between Purdue University and Goddard Space Flight 
Center seeks to further investigate potential links 
between operational libration point orbit stationkeeping 
and any orbit stability information available from dyna- 
mical systems theory. Analyzing the eigenstructure of a 
libration point orbit, a process denoted here as mode 
analysis, supplies information regarding the orbit stability 
as described in Section 3.2. Note that while the mono- 
dromy matrix, M, as discussed previously, is defined 
specifically for periodic orbits, it can be approximated 
for quasi-periodic trajectories, such as those employed for 
the ARTEMIS spacecraft, by integrating the orbit for 
approximately one revolution. Application of mode ana- 
lysis to the ARTEMIS libration point orbits allows for 


a b 
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Fig. 12. ARTEMIS PI EM L 2 stable and unstable 


comparison between the theoretical predictions and 
implemented ARTEMIS stationkeeping maneuver direc- 
tions computed via the OCS during mission operations. 

Using operational ARTEMIS orbit determination solutions 
along with the operational stationkeeping maneuvers that 
were actually implemented (designed using the Optimal 
Continuation Strategy), an approximate monodromy matrix 
is computed by generating and propagating the state transi- 
tion matrix, <£(t 2l Ij), from an initial state for approximately 
one orbital period. The STM is computed in this analysis by 
first adding a position or velocity perturbation - that is, 
4 x 10“ 4 km in position and 1 x 10 -4 cm/s in velocity - to 
each component of the initial state and integrating the 
trajectory forward in time. Then, a finite-difference STM 
using initial and final state information from numerical 
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x (x 10 4 km) 

directions, (a) x-y projection and (b) x-z projection. 



Fig. 13. ARTEMIS PI and P2 EM L, stable and unstable mode directions, (a) ARTEMIS PI and (b) ARTEMIS P2. 
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integration yields an approximation of the monodromy respectively. The stable and unstable eigenvector direc- 

matrix The mode information is generated, Le., the local tions are also displayed for one revolution of the PI EM L, 

stable and unstable eigenvalues and eigenvectors are com- orbit in Fig. 13(a) and for the P2 EM L, orit in Fig. 13(b). 

puted, and the actual maneuver direction is compared with In Figs. 14 and 15. the stationkeeping AV directions and 

the stable/unstable eigenvector information. the stable/unstable mode directions corresponding to 

As an example, the stable and unstable mode direc- the location of a particular maneuver are plotted for the 

tions associated with the ARTEMIS PI EM l 2 orbit over ARTEMIS PI and P2 spacecraft in EM I, orbit. The angle 

one revolution are represented in Fig. 12 in blue and red, between the AV vector, as expressed in the EM rotating 



x (x 10 4 km) 


Fig. 14. ARTEMIS PI SK Maneuver #21. 



x (x 10 4 km) 


Fig. 15. ARTEMIS P2 SK Maneuver #4. 
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Fig. 16. Total (blue) and in-plane (red) angles between AV vector and the associated stable mode direction, (a) ARTEMIS PI and (b) ARTEMIS P2. 
(For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.) 


coordinate frame, and the stable mode direction is com- 
puted from the dot product of the two vectors and is 
presented for all PI and P2 stationkeeping maneuvers in 
Fig. 16. The results of this investigation demonstrate that, 
even when spacecraft maneuver direction constraints are 
considered [1,2], the ARTEMIS stationkeeping maneuvers 
that were executed during mission operations using the 
Optimal Continuation Strategy are generally aligned with 
the stable eigenvector direction. This is an interesting 
result given that previously developed stationkeeping 
methcds incorporating Floquet analysis typically employ 
stationkeeping maneuvers that align with the unstable 
mode in such a way that the unstable component of the 
error is simply canceled. Perhaps not surprisingly, the 
formulation of the OCS apparently yields maneuvers that 
strive to return the vehicle to the vicinity of the path 
dictated by the original libration point orbit insertion 
condition as the most effective approach. Given the 
ARTEMIS post-mission mode analysis, it appears likely 
that stable mode direction information can be incorpo- 
rated into existing stationkeeping algorithms to compute 
optimal libration point orbit maintenance maneuvers 
more quickly and efficiently as demonstrated recently 
by Pavlak and Howell [29]. At a minimum, the stable 
mode direction may serve as an initial guess in an optimal 
strategy. While the full ARTEMIS post-mission mode 
analysis is still being examined, a basic conclusion 
is that maneuver placement along the stable mode direc- 
tion can be used to maintain a libration point orbit by 
effectively perturbing the trajectory back into the center 
subspace. 

6. Summary and conclusions 

The ARTEMIS mission utilized high-fidelity modeling and 
an optimal stationkeeping strategy to successfully maintain 
the PI and P2 spacecraft in Earth-Moon Li and L 2 libration 
point orbits for eleven months. This research seeks both to 
summarize the results of ARTEMIS mission operations and 
to leverage the circular restricted three-body problem and 
dynamical systems theory in an effort to gain further insight 
into the evolution of the ARTEMIS libration point orbits and 
the stationkeeping maneuver design process. The ARTEMIS 
libration point orbit insertion conditions were dictated by 


the ballistic Sun-Earth to Earth-Moon transfers that were 
utilized and, consequently, the ARTEMIS libration point 
orbits were selected in to enable low-cost transfers between 
EM L 2 and Li and, ultimately, into low inclination lunar 
orbit. Investigating the long-term behavior of the ARTEMIS 
PI and P2 orbits using Poincare maps in the circular 
restricted three-body problem reveals that each of the 
ARTEMIS EM libration point orbits are, in fact, arcs of large 
quasi-halo orbits. Analysis of the large quasi-halo orbits 
available in the CR3B problem demonstrates that these 
orbits include both high and low z-amplitude phases and 
are well-suited for the ARTEMIS mission because they 
facilitate both highly out-of-plane and nearly planar arrival 
and/or departure conditions. With an understanding of the 
evolution of the trajectories represented on the Poincare 
map, maps could serve as a useful guide for the selection of 
orbits with specified characteristics in future mission design 
applications. 

While there are a number of stationkeeping strategies 
available that incorporate multi-body dynamics, the actual 
mission applications and constraints must also be consid- 
ered. The Optimal Continuation Strategy is a general sta- 
tionkeeping algorithm that is capable of computing low-cost 
orbit maintenance maneuvers without a reference trajectory 
and is well-suited for maintaining the ARTEMIS Earth-Moon 
libration point orbits while satisfying mission constraints. 
Targeting goals formulated in terms of x-velocity at succes- 
sive x-z plane crossings are utilized to maintain the EM L, 
and L 2 for 1-2 revolutions downstream. Incorporating opti- 
mization, the required stationkeeping AV can be minimized 
and has been demonstrated to be only ~ 5-7 m/s per year 
for the ARTEMIS spacecraft. The optimal stationkeeping 
maneuver directions yielded by the OCS are further ana- 
lyzed by completing mode analysis for the ARTEMIS libra- 
tion point orbits, revealing the directions of the stable and 
unstable eigenvectors. Optimized maneuver directions com- 
puted via the OCS compare closely with the stable mode 
direction for all stationkeeping maneuvers executed, even 
when constraints on the feasible spacecraft maneuver 
directions are considered. With the Earth-Moon libration 
point orbiting phase of the ARTEMIS mission completed, 
investigation continues on robust strategies to improve 
stationkeeping maneuver design for future libration point 
missions. 
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